Extended Ginzburg-Landau formalism: systematic expansion in small deviation from 

the critical temperature 
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Based on the Gor'kov formalism for a clean s-wave superconductor, we develop an extended 
version of the single-band Ginzburg-Landau (GL) theory by means of a systematic expansion in the 
deviation from the critical temperature T c , i.e., r = 1 — T/T c . We calculate different contributions to 
the order parameter and the magnetic field: the leading contributions (oc r 1 ' 2 in the order parameter 
and oc r in the magnetic field) are controlled by the standard Ginzburg-Landau (GL) theory, while 
the next-to-leading terms (oc r 3 ^ 2 in the gap and oc t 2 in the magnetic field) constitute the extended 
GL (EGL) approach. We derive the free-energy functional for the extended formalism and the 
corresponding expression for the current density. To illustrate the usefulness of our formalism, we 
calculate, in a semi-analytical form, the temperature-dependent correction to the GL parameter 
at which the surface energy becomes zero, and analytically, the temperature dependence of the 
thermodynamic critical field. We demonstrate that the EGL formalism is not just a mathematical 
extension to the theory - variations of both the gap and the thermodynamic critical field with 
temperature calculated within the EGL theory are found in very good agreement with the full BCS 
results down to low temperatures, which dramatically improves the applicability of the formalism 
compared to its standard predecessor. 

PACS numbers: 74.20. De, 74.20.Fg, 74.25. Bt, 74.25.Ha 



I. INTRODUCTION 

In 1950 Ginzburg and Landau proposed a phenomcno- 
logical theory of superconductivity (the GL theory) based 
on a specific form of the free energy functional con- 
structed in the vicinity of the critical temperature from 
Landau's theory of second-order transitions^ Minimiza- 
tion of this functional generates the system of two GL 
equations that give the spatial distribution of the super- 
conducting order parameter (the condensate wave func- 
tion) and of the magnetic field in a superconductor. Over 
the years, the GL approach has been enormously success- 
ful in describing various properties of superconductors 
(see e.g. Ref. y), and has been extensively used partic- 
ularly in the last decade in the domain of mesoscopic 
superconductivity, see, e.g., Ref. H 

In 1957 Bardeen, Cooper and Schrieffer (BCS) de- 
veloped the well-known microscopic theory of the con- 
ventional superconductivity- and then, two years later, 
Gor'kov showed that the GL equations can be obtained 
from the BCS formalism^ However, in spite of the avail- 
ability of the microscopic theory, the GL approach re- 
mains still an optimal choice for many practical calcula- 
tions when the spatial distribution of the pair condensate 
and, thus, of the magnetic field are nontrivial, e.g., for 
multiple-vortex configurations. The simple differential 
structure of the local GL equations in comparison with 
the microscopic theory offers clear advantages, including 
the possibility of analytical derivations in many impor- 
tant cases. 

A desire to develop an extension to the GL theory, 
with the idea of improving the formalism while retain- 
ing at least some advantages of its original formulation, 



stimulated significant efforts by many theorists. Several 
GL-like theories of different complexity were proposed. 
In the earliest developments^ the so-called "local su- 
perconductor" formalism was attempted, being a compli- 
cated synthesis of the BCS and GL approaches. The GL 
theory with nonlocal corrections (i.e., including higher 
powers of the gradients of the order parameter) was used 
in studies of the anisotropy of the upper critical field 
(see Ref. d) and the vortex structure (see Ref. @) in d- 
wave superconductors. Recently, various extensions to 
the GL theory were introduced in the context of study- 
ing the Fulde-Ferrell-Larkin-Ovchinnikov state.— ~— In 
all the above examples, extending the GL theory was 
based on the expansion of the self-consistent gap equa- 
tion by including higher powers of the order parameter 
and its spatial gradients phenomenologicallyJ^ However, 
accounting for such higher-power terms is not as straight- 
forward as it may seem. The fundamental problem here 
is to select properly all relevant contributions of the same 
order of magnitude (accuracy). This (serious) issue does 
not arise in the derivation of the original GL theory for 
a single-band superconductor, where only the first non- 
linear term and the second-order (leading) gradient of 
the condensate wave function are included. However, as 
recently shown in Ref. [bH . a similar selection performed 
for the GL theory of a two-band superconductor leads to 
the appearance of incomplete higher-order contributions. 
Such incomplete terms may cause misleading conclusions 
and should be avoided. 

To tackle this problem, one needs to work with a single 
small parameter in the expansion. In the present case, 
such a small parameter is the proximity to the critical 
temperature, i.e., r = 1 — T /T c . Indeed, the standard 
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GL approach can be seen as the lowest-order theory in 
the r-expansion of the self-consistent gap equation, see, 
e.g., Refs. [3 and [H. However, next orders in r are 
also of great importance, for example, to capture the 
physics of different healing lengths of different conden- 
sates in multi-band superconductors . 15 ' 16 In the present 
paper, we show that next orders in r are also important in 
the single-band case, surprisingly improving the GL the- 
ory. We obtain a systematic series expansion of the self- 
consistent gap equation for a single-band, s-wave, clean 
superconductor, using r as the governing small parame- 
ter. In the derivation we employ a technique in the spirit 
of the asymptotic expansion methods used extensively in 
the applied mathematics. 17 Similarly to the asymptotic 
analysis in other models, we obtain a hierarchical sys- 
tem of the so-called transport equations which need to 
be solved recursively starting from the lowest order. Us- 
ing this method, corrections to the standard GL theory 
can in principle be calculated with an arbitrary accu- 
racy. However, unlike asymptotic expansions for linear 
models, here the complexity of the higher-order equations 
increases rapidly and their solution cannot be obtained 
in the general form. Based on the Gor'kov Green func- 
tion formalism, we derive and investigate the first three 
orders of the r-expansion of the gap equation, i.e., r"/ 2 
with n = 1,2,3. To the order r 1 / 2 , we find the equa- 
tion for the critical temperature. Collecting the terms 
proportional to t 3 / 2 , we obtain the standard GL theory 
giving the lowest-order (in r) contributions to the super- 
conducting condensate, i.e., oc t 1 / 2 , and to the magnetic 
field, i.e., oc r. Then, by matching the terms of the order 
t 5 / 2 , we derive equations for the next-to- leading correc- 
tions to the superconducting order parameter and mag- 
netic field, cx t 3 / 2 and r 2 , respectively. The equations 
controlling the order parameter and the magnetic field 
up to the next-to-leading order in t constitute the ex- 
tended GL formalism (EGL). To illustrate the power of 
our extension to the GL theory, we investigate the energy 
associated with a surface between the superconducting 
(S) and normal (N) phases. In particular, we calculate 
the temperature dependence of k* , the value of the GL 
parameter at which the surface energy becomes zero. It is 
important to stress here that contrary to other available 
extensions of the GL theory discussed above, our formal- 
ism is not much more complicated than the standard GL 
theory. As is shown in the calculation of the S /N surface 
energy, plenty of important information can be obtained 
from the EGL formalism analytically. 

Note that as is known, the GL theory is heavily used 
also in the study of thermal fluctuations. Here we do 
not address this issue but deal with an extension to the 
GL formalism in the mean-field level. Our aim is to ex- 
pand the validity domain of the GL theory down to lower 
temperatures, which will be useful for the problems with 
a strongly nonuniform distribution of the pair conden- 
sate, e.g., for multiple vortex solution in the presence of 
stripes and clusters of vortices. In this case any fully 
microscopic approach will be an very complicated and 



rather time consuming task. 

The paper is organized as follows. In Sec. |TT]we intro- 
duce a general approach to construct the asymptotic ex- 
pansion in r for the self-consistent gap equation. In order 
to illustrate the main ideas behind the method, a simpler 
case of zero magnetic field is considered here. Sec. IIIII 
presents the generalization to a nonzero magnetic field. 
Such a generalization requires the normal-metal Green 
function beyond the traditional Peierls (phase) approx- 
imation, and the corresponding expression is presented 
and discussed. The series expansion of the free-energy 
functional up to the next-to-leading order in r is given 
in Sec. IIVI The complete set of equations for the order 
parameter and the magnetic field in the EGL approach 
is derived by finding the stationary point of the func- 
tional in Sec. [V] In Sec. I VII we estimate the accuracy of 
the EGL formalism by comparing its results for the uni- 
form order parameter and the critical magnetic field with 
those of the standard GL approach and the BCS theory. 
Here we demonstrate that the temperature in which the 
GL theory is valid, is dramatically increased due to the 
extension. In Sec. IVIII we investigate the S/N surface 
energy and calculate the temperature-dependent correc- 
tion to k* . Finally, Sec. IVIIII presents the summary of 
our results and our conclusions. 

The main text of the paper contains only the key for- 
mulas, the general ideas and the main steps of the deriva- 
tion. Readers interested in details are referred to Appen- 
dices. In particular, Appendix |A1 shows how to calculate 
the coefficients appearing in the r-expansion of the gap 
equation in the absence of a magnetic field. Appendix [Bl 
presents details of our calculations for the normal-metal 
Green function beyond the Peierls phase approximation. 
Appendix [C] generalizes the calculations given in Ap- 
pendix [A] to the case of a nonzero magnetic field. 



II. SERIES EXPANSION IN r OF THE GAP 
EQUATION AT ZERO MAGNETIC FIELD 

As is known since the classical work by Gor'kov ) 2 i 5 i 18 
the GL equations can be derived from the microscopic 
BCS theory in the most elegant way via the Green func- 
tion formalism. For the sake of clarity of presentation of 
our main ideas, the case of zero magnetic field is consid- 
ered first, while the generalization to a nonzero magnetic 
field is given in the next section. The goal of our work 
is to construct the extended GL formalism through the 
expansion of the Gor'kov equations in r = 1 — T/T c , 
with T c the critical temperature (T < T c ). We start by 
writing the Gor'kov equations as the Dyson equation for 
the Green function in the Gor'kov-Nambu 2— matrix 
representation, which reads as (see, e.g., Ref. Il8l - [20l) 

^ = ^ 0) + ^ 0) A^, (1) 
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with 



given by 



Qu = 



fg u r u \ d(0) ( gL 0) o \ 



where Jvu) = 7rT(2n + 1) is the fermionic Matsubara fre- 
quency (n is an integer and ks is set to unity) and the 
2x2 matrix gap operator A in Eq. ((T|) is defined by 



A = I £ )• (r|A|r')=rf(r-r')A(r'l. (3) 



where the superconducting order parameter A(r) obeys 
the self-consistent gap equation, i.e., 



A(r)=- ff T^JL(r,r), 



(4) 



with g the (Gor'kov) coupling constant. As usual, the 
sum in the r.h.s. of Eq. (j4]) is assumed to be properly 
restricted to avoid the ultraviolet divergence. Equations 
(dJ and © further give 

JUr,r')= |d 3 2/ ^°)(r,y)A(y)^(y,r'), (5a) 
^(r,r') = s£ 0) (r,r') 

+ /d 3 y^°)(r,y) A*(y)^(y,r'). (5b) 



For the normal-state temperature Green function 
G^ (r, y) we have (at zero magnetic field) 



^ 0) (r,y) 



d 3 k e lk ( r ~ y ) 



(27r) 3 ihuj — £k 



(6) 



with the single-particle energy = h 2 k 2 /2m — p, mea- 
sured from the chemical potential /i, and u (r, y) = 

-G- U {y,t). 

The Gor'kov equations ([3]) supplemented by Eq. © 
make it possible to express the anomalous (Gor'kov) 

Green function J- U (r, r') in terms of A(r) and G^ (r, y). 
Then, inserting this expression into Eq. one obtains 
the self-consistent gap equation. Solution to the gap 
equation can be represented in the form of a perturbation 
series over powers of A(r) (which is small in the vicinity 
of the critical temperature T c ): 

A(r)= /d 3 y^ Q (r,y)A(y) + ff[ d 3 % K b (r, {y} 3 ) 

5 

x A( yi )A*(y 2 )A(y 3 ) + K ^{yh) 



3=1 



x A( yi )A*(y 2 )A(y3)A*(y 4 )A(y 5 ) + . . . , (7) 
where {y}„ = {yi, . . . ,y n } and the integral kernels are 



K a (r,y) = -gT £ 0<°>(r, y)Q™ (y, r), 

K b {v,{yh) = -gT ]T ^ 0) (r,yi)^ 0) (yi,y 2 ) 

a; 

x^ 0) (y2,y 3 )^ 0) (y3,r), (8) 

K c (r, {y} 5 ) = -gT ]T <?£ 0) (r, yx)^ ) ( Yl , y 2 ) 

UJ 

x ^ 0) (y 2 , y 3 )^ 0) (y 3 , y 4 )^ 0) (y 4 , y 5 )^ 0) (y 6 , r). 

Equation (J7J can be truncated to a desired order, which 
yields a non-linear integral equation. The latter is further 
converted into a non-linear partial differential equation 
by using the gradient expansion 

A(y,) = A(r + zj) = £ - {z j ■ V r )"A(r), (9) 

n=0 

with Zj • V r the scalar product. In particular, the GL 
equation is obtained when keeping only the first two 
terms, including K a and Kb, in Eq. ([7]) and the second- 
order spatial derivatives in the gradient expansion ([9J. 

While Eq. ([7]) is a regular series expansion of the gap 
equation the partial differential equation mentioned 
above is not. The gradient expansion introduces a second 
small parameter together with the corresponding trunca- 
tion approximation, and the relation between the order 
parameter and its gradients is not known a priory. As a 
result one cannot compare the accuracy of the relevant 
terms in the expansion and the truncation procedure be- 
comes ill-defined. This problem does not appear in the 
derivation of the GL equation where one keeps only the 
lowest second order gradient term. In order to extend 
the GL formalism, one has to deal with a single small 
parameter in the system that can be used to compare 
the relevant contributions. This small parameter follows 
from the solution of the GL equation. When T — » T c 
the order parameter decays as A a t 1 / 2 — > 0. Also the 
solution reveals the scaling length £ oc t -1 / 2 , i.e., the 
GL coherence length, which determines the spatial vari- 
ations of the order parameter in the vicinity of T c and 
dictates that VA cx r, or, using the short-hand nota- 
tion, V oc t 1 / 2 . Thus, the parameter r controls both 
relevant quantities and can be used to produce a single- 
small-parameter series expansion of the gap equation Q. 

A systematic expansion of the gap equation in r can 
be facilitated by introducing the scaling transformation 
for the order parameter, the coordinates and the spatial 
derivatives of the order parameter in the following form: 



J/2 A, 



(10) 



Note that in terms of the scaled coordinates, the typical 
spatial variation of the order parameter occurs on a scale 
that is independent of r to the leading order. After the 
transformation given by Eq. (|10[) . the parameter r enters 
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the expansion in Eqs. (J7J and (J9)) as follows. In Eq. ([9]), 
only coordinate r changes. The scaling of the difference 
z does not change the expressions as it is an integration 
variable, and thus the scaling will not be applied to it. 
As A is now a function of r, each derivative V in the 
expansion ([9]) adds a factor r 1 / 2 . Inserting Eq. (|9]) with 
these factors into the expansion given by Eq. ([7J and tak- 
ing into account the factor r 1 / 2 for the order parameter 
in Eq. (jlOp we arrive at a simple mnemonic rule to count 
the minimal order of each term in the expansion of the 
gap equation: each occurrence of A or V in the formulas 
adds the factor t 1 / 2 . The final form of the expansion is 
obtained by calculating the relevant coefficients through 
the evaluation of the remaining integrals. As those coeffi- 
cients depend on temperature, they can also represented 
as series in r. The formulated procedure allows one to 
calculate the r-expansion for the gap equation to arbi- 
trary order. However, in practice, calculations of higher 
orders become more and more complicated. In this work 
we limit ourselves to the analysis of the self-consistent 
gap equation in the first three orders, i.e., up to the order 
t 5 / 2 . Collecting terms of the order r 1 / 2 , we obtain the 
equation for T c . Working in the order r 3 / 2 , we recover 
the standard GL theory producing the leading contribu- 
tion to A, i.e., ex t 1 / 2 . The order t 5 / 2 yields the equation 
that controls the next-to- leading contribution to A, i.e., 
oc t 3 / 2 (this is what we call the EGL formalism) . Details 
of the selection of all the necessary terms in Eq. ([7]) that 
contribute to one of the three orders mentioned above, 
are given in AppendixlA"! The final result reads (V = Vf ) 



Euler constant and CO ■ •) is the Riemann zeta-function. 
It is of importance to note that Eq. (fTTj) contains only 
half-integer powers of r. The reason for this is two-fold. 
First, due to the structure of Eq. ([7]), there appear only 
odd integer powers of the order parameter. Second, the 
spherical symmetry dictates that the integrals with an 
odd number of V operators are equal to zero. 

The solution to the gap equation must also be 
sought in the form of a series expansion in r. Based on 
Eqs. (Tnj) and (fT2|) . we are able to introduce 



A(r)=A (r) + rA 1 (r) 



(13) 



Substituting this into Eq. (|11[) and collecting terms of the 
same order we obtain a set of equations for each A n . 
Collecting terms of the order r 1 / 2 , we obtain 



-A T )A = 0. 



(14) 



The solution to this equation, i.e., qAt = 1 gives the 
ordinary BCS expression for the critical temperature, i.e., 
T c = (2e^/n)huj D exp[-l/ 5 A(0)]. 

In the order r 3//2 one recovers the standard GL equa- 
tion for the leading contribution to the order parameter 
A : 



aA + o|A | 2 A - /CV 2 A = 0. 



(15) 



Note that its standard form with the temperature depen- 
dent a-coefhcient is obtained by multiplying all terms by 
the factor t 3 / 2 and returning to the unsealed quantities, 



.1/2 



A = air 1 / 2 A + a 2 r 3 / 2 V 2 A + a 3 r 5/2 V 2 (V 2 A) 



-6ir 3 / 2 |A| 2 A - 6 2 t 5 / 2 |2A |VA| 2 + 3A*(VA) 2 

+A 2 v 2 A* + 4|A| 2 v 2 A] + cir 5/2 |A| 4 A, (11) 

where the coefficients dj, b-i and Ci are obtained from the 
integrals with the kernels K ai Kf, and K c , respectively, 
and they are given by 



ai = ^T-a(r +T + 0(^)) > _ = In 
6 1 =6(l + 2r + 0(r 2 )), b=N(0)^)_, 

*=c(l + (r)), 0=^(0)^, 

a 2 =/C(1 + 2t + 0(t 2 )) , K = \h 2 v 2 F , 

6 

03 = C(l + 0(r)), Q=±tfv%, 
b 2 = £(l + (D(T)), £=^h 2 v 2 F , 



(12) 



where a = -N(Q) and N(0) = mk F /(2ir 2 h 2 ) is the DOS 
at the Fermi energy, with vp the Fermi velocity; ojp) de- 
notes the Debye (cut-off) frequency, 7 = 0.577 is the 



arAo + 6|A | 2 A - /CV 2 A = 0. 

Finally, collecting all terms of the order t 5 / 2 , we arrive 
at the equation for Ai, i.e., the next-to-leading term in 
the order parameter: 

0A1 + 6(2|A | 2 A! + A 2 Al) - KV 2 A l = F, (16) 

where F is given by 

F = - |Aq + 2/CV 2 A + QV 2 (V 2 A ) 

- 26|A | 2 Ao - J2A0 I VAo| 2 + 3Aq (VA ) 2 



A 2 V 2 A* + 4|A | 2 V 2 A C 



c|Anl 4 An 



(17) 



This is a linear differential inhomogeneous equation to be 
solved after Ao is found from Eq. fp~5|) . Note that similar 
features in the r-expansion of the gap equation appear 
for a two-band superconductor, as welli 1 ^ 

We note again that in principle, one can continue the 
procedure up to arbitrary order in r, obtaining correc- 
tions to the standard GL theory with desired accuracy. 
While the equation for Ao is nonlinear, the higher order 
contributions to A will be controlled by inhomogeneous 
linear differential equations. Such a system of equations 
is solved recursively, starting from the lowest order, since 



5 



solutions for previous orders will appear in the higher or- 
der equations, but not vice versa. The solution to the 
system will thus be uniquely defined (when the relevant 
boundary conditions are specified), ensuring consistency 
of the developed expansion. 

We also remark that the structure of Eq. (p~6|) makes it 
possible to conclude that the next-to- leading term Ai(r) 
is not trivially proportional to Ao(r). For that reason, 
the spatial profile of Ai(r) is different compared to A (r). 
This means that the characteristic length for the spatial 
variations of A(r) in EGL differs from the standard GL 
coherence length. However, both lengths have the same 
dependence on r, i.e., VAi oc r. 



III. SERIES EXPANSION IN r OF THE GAP 
EQUATION FOR NONZERO MAGNETIC FIELD 

The magnetic field enters the formalism developed in 
the previous section via changes in the normal-metal 
Green function. In the Gor'kov derivation the field- 
induced corrections are accounted for through the field- 
dependent Pcicrls phase factor as 



^ r (rt,rY) = e 



iM-dq 



(18) 



where A is the vector potential, i.e., B = rot A. It is 
of importance to note that the integral in the exponent 
is evaluated along the straight line connecting r' and r. 
This approximation leads to Eq. (fT5|) . where the gauge- 
invariant derivative replaces V. Obtaining corrections to 
the GL equation requires the Green function to be cal- 
culated with an accuracy sufficient to produce the com- 
plete set of terms contributing up to the order r 5 / 2 in 
the r-expansion of the gap equation. Taking into ac- 
count Eqs. CL!} and H2|), one concludes that the normal- 
metal Green function must be calculated with the accu- 
racy C(r 2 ). 

Accounting for the magnetic field in the expansion can 
be done by noting that the critical magnetic field H c in a 
superconducting system is proportional to r. Similarly a 
solution of the GL equations for the field also changes as 
oc t. Thus the derivation of the r expansion for the field 
corrections can be conveniently done by introducing the 
following scaling for the magnetic field: 



A = r 1/2 A, B = rotA = tB, 



(19) 



so that the critical field quantities become independent 
on r in the first order. We note that the spatial depen- 
dence of the magnetic field can also be represented in 
the form of the gradient expansion as it is done for the 
order parameter in Eq. ((9]). As the characteristic scale 
for the spatial variations of the magnetic field is defined 
by the magnetic penetration depth A cx r -1 / 2 , which has 
the same r dependence as the GL coherence length £, we 
can again employ the scaled spatial coordinates as those 
in Eq. (| lOj) . So, the gradient expansion for the magnetic 



field follows the same rule as for the order parameter: 
each V introduces an additional factor t 1 / 2 . 

The r-expansion of the Green function can now be 
done using the path integral method, by accounting of 
the quantum fluctuations around the classical trajectory. 
For details of the calculations we refer the reader to Ap- 
pendix [Bl The final expression for the Green function 
that contains all contributions to the order r 2 reads as 



CTr,r')=e 



/A-dq 



B 2 (r) 



24m 2 c 

5/2^ 



I + l m (r-r')ld u \ +0(r 5 / 2 )| G<% =0 (r,r'), 

(20) 

where the phase factor in the exponent is the same as in 
Eq. (|18l) . We remark that the Peierls phase also contains 
terms of higher orders than t 2 . However, they do not 
enter the final equations for the next-to-leading contri- 
bution to the order parameter and magnetic field. It is 
simply convenient to represent the Green function in the 
form with the Peierls factor because it naturally leads to 
the appearance of the gauge invariant spatial derivatives 
of the order parameter. One also notes that this factor is 
written using the unsealed quantities. In fact, the scaling 
does not affect it. In addition, we do not scale the differ- 
ence of the coordinates z = r — r' as it is an integration 
variable. The new field-dependent term in Eq. (|20l) is 



proportional to r B . It does not follow from the Peierls 
phase and, so, is not present in the approximation given 
by (fl"8]) . It is interesting that the field gradients do not 
contribute to the correction to the Peierls approximation. 

The field-modified expansion of the gap in powers of 
the order parameter and its spatial derivatives are ob- 
tained by substituting Eq. ([20]) into the kernels in Eq. ([8|) 
and proceeding in a manner similar to that discussed in 
Appendix lAl for Eq. (fTTj) . It is convenient to remove the 
phase factor from the Green functions and introduce the 
"two-point" auxiliary order parameter, i.e., 



A(f,f') 



-ffM-dq _ 
e f A(f). 



(21) 



Then, as shown in Appendix \C\ the modification of 
Eq. (fTT|) due to the phase factor is that A(r) is replaced 
by A(r,r'), and the limit r' — s- r is implied after all rel- 
evant calculations (we remark that such a limit is not 
permutable with the differentiating with respect to r). 
In particular, for the first three terms in the modified 
Eq. (fTTj) we have 



a^ 1 ' 2 lim A(f,f') = air 1 / 2 A, 

r'— ¥v 

a 2 r 3/2 lim V 2 A(r,f') = a 2 T 3/2 D 2 A, 

r' — »r 

a 3 r 5/2 lim V 2 (V 2 A(r, ?')) = 



4ie _ .- - 4e 2 



(22a) 
(22b) 

A, (22c) 
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with D = V — iff A, the gauge invariant gradient, and 
rotB = V x B. The terms related to the kernel K b in 
the field-modified Eq. (fTT|) read 



-6ir 3/2 lim |A(f,r')| 2 A(r,f') = -6it 3/2 |A| 2 A, (23a) 



theory^ Its expansion in A reads 



-6 2 r 5/2 lim 



2A(r,f') |V r A(r,f')| 2 +3A*(r,f') 



(V r A(r,f')) +A 2 (f,f') V 2 A*(f,f' 



-4|A(r,r')| 2 V 2 A(r,r') 



= -fe 2 r 5 / 2 



2A IDAI 



- 3A*(DA) 2 + A (D A)* + 4|A| 2 D 2 A , (23b) 



whereas the term coming from the integral with the ker- 
nel K c is of the form 



cir 5 / 2 lim |A(r,r')| 4 A(r,f') =cir 5 / 2 |A| 4 A. (24) 



In addition to the contributions that appear due to the 
Peierls factor in the Green function, we also obtain an 
extra contribution to the r.h.s. of Eq. (fTTj) due to the 
terms proportional to B 2 in Eq. (|2T)]) . This contribu- 
tion comes only from the integral involving the kernel 
K a (when keeping terms up to the order t 5 / 2 ) and reads 



-a 4 r 5 / 2 B 2 (f) lim A(f , r') = -a 4 r 5 / 2 B 2 A, (25) 



where a 4 = &fi, 2 e 2 /(36m 2 c 2 )(l + 0(t)), with b given by 

Eq. CE2D. 



Using the modified Eq. (fTTj) together with Eqs. 
(1241) and collecting terms of the same order in t, one can 
generalize Eqs. (|T5j) and (fT6|) to the case of a nonzero 
magnetic field. However, since the magnetic field is also 
a variable in the GL theory, it needs to be found sclf- 
consistently from a complementary set of equations. The 
most elegant way to derive the complete set of equa- 
tions for the magnetic field and for the order parameter is 
based on the free-energy functional that accounts for the 
energy associated with the presence of the magnetic field 
and the superconducting pairing. This functional must 
be also represented as a series expansion in r, which is 
the subject of the next section. 



IV. FREE-ENERGY FUNCTIONAL 



The free-energy functional J- s can be obtained, e.g., by 
using the path integral methods developed for the BCS 



T, = F n ,B=o + / d 3 r 



B 2 (r) 

8tt 



i Jd 3 rd 3 y [S(r - y) - K a (r,y)] A*(r)A(y) 

3 

^ |d 3 r [] d 3 ^ K b (r, {y} 3 ) A*(r)A( yi ) 
A*(y 2 )A(y 3 ) - i jd'-r\l\ K c (r,{y} 5 ) 



x A*(r)A( yi )A*(y 2 )A(y 3 )A*(y 4 )A(y 5 ) 



(26) 



where T n denotes the free energy of the normal state 
(-7 r n,s=o stays for the zero magnetic field). Here it is 
worth noting that, generally, 

K* a (r,y) = K a (y,r), 

#b(c)( r >{y}3(5)) = ^6(c)({y}3(5),r), 

which means that the r.h.s. of Eq. (|26l) is a real quantity 
(as necessary for the free energy). In addition, we have 

K b (r,yi,y 2) y 3 ) = K b (y 2) y 3 , r, yi), 
^ c (r, y i,y 2 ,y 3 ,y 4 ,y 5 ) = ^ c (y2,y3,r,yi,y4,y5) 
= ^ c (y2,y3,y4,y5,r,yi), 

which makes it possible to immediately find that the ex- 
ternum condition of this functional with respect to A* 
leads to Eq. 0. 

Series expansion of the free-energy functional given by 
Eq. ([2"fJ| is obtained by following essentially a similar ap- 
proach as is used in previous sections. The scaling trans- 
formation for the order parameter, coordinates and mag- 
netic field is introduced, and the gradient expansion for 
the order parameter © is substituted into Eq. (|26|) . As 
before, the coefficients of the series expansion appears as 
the integrals with kernels K a ^ >c . As we are interested 
in the terms that are used to derive the GL equations 
and the equations for the next-to-leading contributions 
to the order parameter and the magnetic field, we must 
expand the functional up to the order t 3 . This follows 
from the rules for counting the powers of r introduced 
in Sec. |TTJ A calculation of the series expansion in r of 
the free-energy functional is greatly simplified by noting 
that all terms in Eq. (|26j) can be derived from the cor- 
responding terms in Eq. (Ill[) by multiplying the latter 
by (2/n)A*, where n is the number of A*-factors of the 
corresponding integral term in Eq. (|26l) . As a result, we 
obtain the functional, which can be represented in the 
following symmetric (real) form (here and below we omit 
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bars over the scaled quantities unless it causes confusion): 
B 2 1 , 



Is — fn,B=0 

- T~a 3 



8tt 

(|D 2 A| 2 ^ 



[g '-ai)|A| 2 +a 2 |DA| 2 

4e 2 



-rotB • i 

3 



h 2 ( 



r B 2 |A| 2 



r a4 B 2 |A| 2 + ||A| 4 -r| 



(A*) 2 (DA) 2 + A 2 (D*A*) 2 



|A| 2 |DA| 2 
-rflAI", 



(27) 



where f s r n \ = .F s ( n )/V (with the scaling / = /r 2 ) and i 
is given by 

i = i|^(AD*A*-A*DA). (28) 
fie 

We stress that i is not the current density j. However, 
when replacing A — > Ao and A — > Ao in Eq. (|28p . we 
find io being proportional to jo, the leading contribution 
to the current density j, see the next section. We also 
note that the representation of the functional in a real 
symmetric form as in Eq. (|27[) implicitly relies on the 
disappearance of the corresponding surface integrals. It 
is easy to check that this is ensured by the boundary 
conditions discussed in Sec. El 

Obtaining the final series expansion in r for the free- 
energy density requires the r-expansion for the coeffi- 
cients given in Eqs. ([12")) and (|25l) . for the order param- 
eter given by Eq. ([13")) , and for the magnetic field. The 
latter is expressed in the form 



A = A + rAi 



B = B +tBi 



(29) 



After collecting the relevant terms, the free-energy den- 
sity f s is written in the form 

fs- f n , B =o = / + r/i+O(r 2 ), 



where the leading-order term (the standard GL func- 
tional) is specified by 



/o = ^ + a |A | 2 + ^|A | 



/C|D A | 



(31a) 



and the next-to-leading contribution /i can be written 



as a sum of two terms, i.e., f\ = f\ 



(0) 



f[ l) with 



fi 



"'-f|A | 2 



2/CIDoAol 



2 A | 2 



1 4e 2 
+ -rotB -i + -^r i B 2 l |Ao| 



fi(|D 

b e 2 h 2 



:77 — t^tBo|Ao| 2 
36 m l t* 



b\A \ 



£ 

2 

*\2 



|Ao| 2 |D Ao| 2 



(AS) 2 (D A 



+ A 2 (D*A*) 



A n | 6 



(31b) 



and 



B Bi 



4tt 



(a + 6|A | 2 )(ASA 1+ AoAt; 



+ K 



(D A -D5A*+c.c.) - Ai -i 



(31c) 



Here D is equal to D with the substitution A — > A , and 
io is obtained from Eq. (|2"51) by A — > A and A — > A . 



V. COMPLETE SET OF EQUATIONS FOR B / 

A system of coupled equations for both A and A are 
obtained from the extremum condition of the free energy 
given by Eq. (|2~T1) . The subsequent expansion of the ob- 
tained expressions yields the GL equations as well as the 
equations for the next-to-leading contributions to the or- 
der parameter and the magnetic field. Alternatively, this 
full set of equations can also be obtained by finding the 
extremum of the free energy with respect to Ao and Ao 
in such a way that all terms appearing in Eq. (|3"0)) are 
taken separately. Note that finding the extremum of the 
functional with respect to Ai and Ai yields the same set 
of equations. This property follows from Eqs. (1131) and 

By searching for the minimum of the functional with 
the density /o, we reproduce the standard GL equations: 



aA + 6|A | 2 A - /CD 2 A - 0, (32a) 



4-7T 

rotBo = — jo, 
c 



(32b) 



where jo = ICciq. Likewise, the same equations are ob- 
tained by finding the extremum of f\ with respect to A\ 
and Ai. 

The minimum of the functional with the density /i 
with respect to Aq and Ao yields the following equations: 

aAi + 6(2 |A | 2 Ai + AgA*) - /CD 2 Ai = F, (33a) 



rotBi = — ji. 
c 



(33b) 



(30) The r.h.s. of Eq. (|33ap is given by 



F 



-A + 2/CD 2 A + Q D 2 (D 2 A 



■ 4e „ ^ , 4e 2 
? 3^c r ' + ft 2 ? 



■B 2 A 



r B 2 Ao-26|Ao| 2 A -£ 



36m 2 c 2~u-u ^-u.-u ~ L 2A |DoAo| 

+ 3A5(D A ) 2 + A 2 (D 2 A )* + 4|A | 2 D 2 1 A 

+ c|Ao| 4 Ao-^/C{A 1 .D } + Ao, (34) 

with {Ai • Do}+ denoting a symmetrized product. The 
next-to-leading contribution to the current density j i [j = 
jo + 7"ji+0(T 2 )] appearing in the r.h.s. of Eq. (|33bl) splits 
into two terms, i.e., 



ji = /Cdi + J, 



(35) 



where 



.2e 



H = i— (A D*At + AiD^A* - AJD A 

8e~ 



A*D A!) -— A^Aol 2 , (36a) 
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and 



J = c{ (2JC - 3£|Ao| 2 )io + Qi' + ^rot rot i 



Q 



h 2 c 2 



rot (B 1 A 1 2 ) - -|A | 2 rotB 



b e 2 h 2 
18 m 2 c 2 



rot(B |A | i 



(36b) 



with 



2e r 



Ao(D D 2 Ao)*-D A (D 2 A )* 
+ D 2 A (D A )*-A*D D 2 A 



(36c) 



We note that, based on Eq. 



one finds that in 



(2//C)(a + 6|A | 2 )i . We also note that ii given by 
Eq. (|36a[) is the next-to-leading contribution to i = 
io + rii + 0(t 2 ). Equations (j3"3"|) - (f3l)|) are the general- 
ization of Eq. (fT6"|) to the case of a nonzero magnetic 
field. 

Through our derivation, we used the straightforward 
generalization of the Ginzburg-Landau boundary condi- 
tions at the specimen surface, i.e., 



D 0± A = 0, B 0± A 1 



. 2c . 
i — Ai_lA 

he 



0. 



(37) 



As seen, Eq. (|3T|) follows from the expansion of Di A = 
in t. These boundary conditions cancel the surface in- 
tegrals appearing in the procedure of the variation of 
the free-energy functional with respect to A and A*. In 
addition, Eq. (|37| allows us to cancel the surface inte- 
grals that appear due to the obvious requirement of the 
self-conjugation of the free-energy functional. Note that 
Eq. (|3T|) is not the only possible choice for the boundary 
conditions that makes all the relevant surface integrals 
equal to zero. While being of great interest, discussion of 
different possible variants of the boundary conditions for 
different interfaces is beyond the scope of the present in- 
vestigation. Below, based on our formulation of the EGL 
formalism, we investigate properties of a bulk supercon- 
ductor that are not sensitive to the particular choice of 
the local boundary conditions at the specimen surface. 



VI. VALIDITY DOMAIN OF THE EGL 
FORMALISM 



In this section we estimate the domain of the quanti- 
tative/qualitative validity of the GL approach when ex- 
tended to the next-to-leading order in r. Obviously, a 
detailed analysis of the accuracy of the EGL formalism, 
including spatially nonuniform solutions of Eq. (f33|). re- 
quires much effort and is beyond the scope of the present 
work. However, it is possible to get a feeling about the 
accuracy of the formalism in question on the basis of 




FIG. 1: (Color online) The temperature dependent gap (un- 
sealed) in units of the zero-temperature order parameter cal- 
culated within the full BCS approach Abcs (0) versus the rel- 
ative temperature T/T c : the solid curve represents the full 
BCS; the dashed curve shows the result of the EGL formalism 
given by Eq. (139[) : the dotted curve illustrates the standard 
GL approach. 



the spatially uniform case. Below we compare the EGL- 
results for the order parameter and the thermodynamic 
critical field H c with those of the BCS model. 
For the spatially uniform case Eq. fp~6|) yields 



Ai 
A 



bulk 



3 
4 



ac 
262 



= -4 (1 - 



31C(5) \ 
49C 2 (3)/ : 



(38) 



with A = {-a/b) 1 ' 2 , the solution of Eq. $To$). Taking 
into account Eq. (fT3")) and using Eq. (|3"5)l . we obtain the 
order parameter in the unsealed representation up to the 
order r 3 / 2 as 



A(T) 
Abcs(0) 



7C(3) 



-1/2 



31C(5) 
49C 2 (3) 



(39) 



where Abcs(0) = (7r/e 7 )T c is the zero-temperature gap 
calculated from the full BCS formalism, see, e.g., Ref . Il8l 
Results found from the standard and extended GL ap- 
proaches are compared to the full BCS solution in Fig. [T] 
We can see that the GL result notably differs from the 
BCS curve below temperatures T = 0.7-0.8T c . At the 
same time the EGL approach is in a very good quantita- 
tive agreement with the BCS theory down to T = 0.2T C , 
and only below this temperature the order parameter cal- 
culated within the extended formalism exhibits a slight 
decrease not supported by the BCS picture. 

The thermodynamic critical field H c measures the con- 
densation energy so that 



TT~ — fn,B=0 — fs,B=0, 



(40) 



where f s .B=a is the free-energy density of a homogeneous 
superconducting state in the absence of a magnetic field. 
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FIG. 2: (Color online) The thermodynamic critical magnetic 
field H C (T) (unsealed) in units of // c ,bcs(0) versus the relative 
temperature T/T c : the solid curve represents the BCS theory; 
the dashed curve shows the result of the EGL formalism; the 
dotted curve illustrates the standard GL theory. 



Using Eqs. (fT5|). (|T5|) . (PU|) and (f?T|). we find 
H c = H c0 + tH c1 +O(t 2 ), 



n _ / 47m2 u 



The solution H c = H c o recovers the result of the GL the- 
ory (H c = tH c q in the original unsealed variables). The 
term tH c \ (t 2 H c i in the original unsealed variables) pro- 
vides the next order correction. The numerical coefficient 
ac/(36 2 ) is calculated using Eq. (TT2"1) . This yields 



= 1 - 0.273t + O(t 2 ). 



(42) 



Being back to the original unsealed variables, we get for 
the thermodynamic critical field up to the order t 2 



H C (T) 



H, 



c,BCS 



(0) 



7£(3) 



31C(5) 



2^ 4!)C-I3'l 



(43) 



where # c , B Cs(0) = ^vr^O)] 1 / 2 7rT c /e 7 is the zero- 
temperature thermodynamic critical field.— Figure [2] 
shows the result given by Eq. (|4"3"|) as compared to those 
of the standard GL formalism and the BCS approach. 
Here one notes a very good quantitative agreement be- 
tween the EGL formalism and the BCS theory down 
to temperatures 0.3-0.4T c . In particular, the EGL re- 
sults are only by 5% larger at T = 0.35T C . At lower 
temperatures the curve representing the EGL formal- 
ism deviates from the BCS data by 10-20%. This can 
be compared with the standard GL theory for which 
HM/Hc,bcs = eV8/(7C(3)) = 1.736. 

Thus, based on the results given in Figs. Q] and [2j 
one is able to expect that the domain of the quantita- 
tive validity of the EGL theory (in the clean limit) is 



r < 0.7 (T/T c > 0.3), which is a significant extension as 
compared to r < 0.2-0.3, typical for the standard GL 
approach. 



VII. SURFACE ENERGY IN THE 
NEXT-TO-LEADING ORDER IN r 

The EGL formalism given by Eq. (|3"3"]l allows one to 
find corrections to the solutions of the physical problems 
for which the GL approach is relevant. Here we illus- 
trate the power of the EGL formalism by investigating, 
in the next-to-leading order in r, the energy associated 
with a surface separating the superconducting and nor- 
mal phases. It is one of the fundamental problems in the 
theory of superconductivity in view of the fact that su- 
perconducting materials are classified as type I or type 
II according to whether the surface energy is positive or 
negative, respectively. From the standard GL theory it 
is well-known that the surface energy is controlled by the 
GL parameter n = A/£, where A is the magnetic pen- 
etration depth and £ is the GL coherence length. The 
surface energy is positive for k < n* and negative for 
k > ft*, where k* = l/v2 is a universal constant, being 
independent of temperature. Now, based on the EGL 
approach, it is interesting to check whether or not k* is 
independent of temperature in the next-to-leading order 
in t. 

Following the standard calculation of the surface en- 
ergy, see, e.g., Ref. l22l we consider a surface separating 
the superconducting and normal phases perpendicular to 
the z-axis. The superconducting and normal phases are 
found at z — > — oo and z — ¥ oo, respectively. This means 
that the superconducting order parameter approaches its 
bulk (uniform) value at z —> — oo and goes to zero at 
2 — > oo. The magnetic field is chosen in the y-direction. 
It approaches the thermodynamic critical field H c [see 
Eq. P4j) ] at z — > oo and becomes zero at z — > — oo. Note 
that z — > ±oo means here that the point is far beyond 
the surface but still inside the specimen. When going 
far beyond the specimen, the magnetic field always ap- 
proaches H c . Further, the vector potential is taken in 
the Coulomb gauge in the form A = {A(z),Q, 0}. The 
magnetic field respectively reads as B = {0, A'(z), 0}. 
Hereafter the prime sign denotes the derivative with re- 
spect to z. It should be noted that this choice applies to 
both Bo (Ao) and Bi (Ai). In this case Eqs. (|32ap and 
(|33ap contain only real coefficients and, hence, have real 
solutions. 

For our choice the surface energy per unit area <j sn 
(the surface tension) is given by a ID integral of the dif- 
ference between the Gibbs free-energy densities of the 
nonuniform superconducting solution g s and the uniform 
normal state g n , i.e., 

f HB 

@ sri = dz (g s - g n ), g s{n) = / s(n) - -j— , (44) 
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where f s is given by Eqs. (|30)) and (j3"Tj) and /„ = 
fn.B=Q + H 2 /(87r). The external magnetic field H is con- 
stant and its absolute value is equal to the thermody- 
namic critical field H c . In the normal phase B = H and 
B = H = H c , which results in g n = f n ,B=o — H 2 /(8n). 
To proceed further, it is convenient to introduce the fol- 
lowing dimensionless quantities: 



A' H c0 \' 
~ _ . 4ttO „ 

lO — lO— ) 9s(n) 



B X b A 

TJ— > A = V A , 

H c0 V a 



H, 



cO 



47T 
n c<A 



9s(n)i &sn — 



47T 



(45) 



<7i can be transformed into the functional derivatives 
of /o [see Eq. pial) ] with respect to A and A . These 
derivatives are equal to zero and generate the standard 
GL equations (|32a[) and (|32b[) that are now reduced to 



A 2 

^A = 0, 



A" 



-io = A A 2 , 



(47a) 
(47b) 



with the boundary conditions (inside the sample) 

A (-oo) = l, A (-oc) = 0, 

Ao(oo) = 0, A'o(oo) = 1. (47c) 



with A = (hc/\e\)y/— b/ (327r/Ca). In what follows, we 
omit the tilde in all the formulas unless it causes confu- 
sion. Using the new dimensionless quantities, we arrive 
at 



a Q + rai, a l 



-t-oo 

d2 9i(z), 



(46a) 



where go(z) and g± (z) arc the coefficients of the expansion 
of the difference g s — g n in powers of r, i.e., g s — g n — 
9o + 9i T + C( t2 )- Taking into account Eqs. piaft and 
(gU), we obtain 



9a 



AqAq 



12 

^0 



A 4 
^0 



1 



1 )Al + ^ + M-l) 



(46b) 



From Eqs. ()31bp . (131c[) and (|44[) . it is seen that g\ can be 
split into two parts, i.e., 



9i — 9\ + 9\ , 



(46c) 



where g^ includes only the quantities with the index 0, 
i.e., 



9i 



(0) 



A 2 
2 

1 

3k 
Ca 



i \2 



2 A^ 

/ \2 



(A' ) A2 



K 2 l\ k 2 2 
(A' ) 2 A 2 



A8(k F \y 



A 4 



H(Ag) 2 A 2 + 3A 2 A 4 ]+^Ag', (46d) 



Here we stress again that far outside the specimen we 
always have A' Q = 1 and A' x = — ~ — However, inside 
the sample, deep in the superconducting domain, we are 



able to employ A' 



0(1) 



(in addition, due to our choice 



of a real order parameter, we also have ^4o(i) = 0)- So, 

based on Eqs. (|46| and (|47|) . we find that g± can be 
replaced by 



9 i" = 



G+£)w-i>- 



(48) 



As seen, there are only two physical parameters that en- 
ters the relevant expressions for <r s „, i.e., the GL param- 
eter k and the product of the Fermi wavenumber and 
the magnetic penetration depth kp\. As to the quanti- 
ties Qa/K. 2 , ac/(36 2 ), Ca/{bJC) and H c i/H c0 , they are 
simply numbers. From Eq. (|12l) we obtain 



Qa 
K 2 



ac 

~°- 817 < 



Ca 

-°- 227 < bJC 



-0.454 



(49) 



and, in addition, H c i/H c q = —0.273, as seen from 
Eqs. (|41~T) and (|42|) . For conventional superconductors 
the term including kp\ in Eq. (|46d[) is extremely small 
and, so, we are left with only one governing physical pa- 
rameter k in both the leading and next-to-leading orders. 

Now, we have everything at our disposal to calculate 
k* at which the surface energy a sn becomes zero, i.e., 



Q ; o(k*) + rai (k*) = 0. 



(50) 



The solution to Eq. (|50p should be represented in the 
form of the r-expansion 



k* +tkI + 0(t 2 ), 



(51) 



and g[^ reads 

g^ =2[(A 2 -l)AoA! 



(i) _oT/a2 ^ a _a. , ^AoAx 



rA Ai 



-A 1 i Q + {A' Q -l)(A' x + \ + ^), (46e) 

incorporating Ai and A±. 

When calculating a sn , we use Ag^) and Ag^) obtained 
by solving Eqs. f[32]) and (|3"3)) . This makes it possible to 
simplify the problem significantly, excluding the terms 
involving Ai and A\. The point is that such terms in 



where Kg obeys 



a {Ko) = 0. 



(52) 



When Kg is known, the next-to-leading order contribution 
k\ is found by expanding Eq. (|50p in powers of r, i.e., 



ao(K.g)+™l (Kg) + tk\ 

h-oo 



dap 



J Z \SA <9k* SA 8k* 



0, (53) 



11 



where all contributions are calculated at k — Kq. The 
partial derivative in the above expression accounts 
for a change in ao caused by the explicit presence of k 
in Eq. (|47al) . The functional derivatives of ao appear in 
the integral in the r.h.s. of Eq. ([53]) due to changes in A 
and Aq when changing n. At first sight, this complicates 
the problem enormously because of the need to know 
the derivatives and at the point n — Kq. It is, 

however, clear that = and = because they 
yield the GL equations in the leading order in r. So, k\ 
is given by 

*i = , (54) 

and the derivatives of A and A with respect to k at the 
point k — k,q disappear from the final result. 

Thus, to calculate both kq and Ki, one only needs to 
find Ao and Aq from the standard GL equations (|47a[) 
and (|47b|) at k = Kq. The calculation of Kq = l/\/2 is 
a classical problem that can be found in textbooks on 
superconductivity, see, e.g., Ref. |22j. Calculating m, we 
first find 

+ 00 
— oo 

where Eq. (I47a[) was used to simplify the expression. In 
turn, using Eqs. (|46|) . (|47|) and the helpful relation 22 

A 2 = 1 - A' (56) 

and assuming kp\ 3> 1, which is always satisfied in the 
conventional superconductors, one finds 

/ 2Qa ac \ 

(2Ca 5Qa ac \ 

+ l2 \-bic~3io~w)' (57) 

with 

+oo 

1 2 = /dzA^l-A 2 ,). (58) 

— oo 

Numerically solving Eqs. (j47|) . we obtain Ii = 0.775 and 
1 2 = 0.433. Finally, substituting Eqs. (JS3J) and (J5TJ) into 
Eq. (JS11) and using Eq. gSJ), we obtain k\ = -0.027^. 
Thus, Eq. (ISTj) for k* reads 

k* = -0.027t + O(t 2 )). (59) 

v2 

As seen, contrary to the result of the standard GL the- 
ory (the leading order in our r-expansion) , k* given by 
Eq. (|59p is temperature-dependent. This means that the 
traditional classification of type-I and type-II supercon- 
ductors becomes, in principle, temperature-dependent. 
However, inconvenience caused by such dependence is 
not crucial because of the very small value of «*, i.e., 
a change in k* with temperature is less than 2% in the 
strict validity domain of the EGL theory. 



VIII. CONCLUSIONS 

Employing an approach similar to the asymptotic- 
expansion methods used in the mathematical physics, we 
constructed a systematic expansion of the self-consistent 
gap equation (for a clean s-wave single-band supercon- 
ductor) in powers of r, the proximity to the critical 
temperature. The procedure of matching the expansion 
terms of the same order of magnitude generates a hi- 
erarchy of equations. The lowest-order theory, i.e., the 
equations for the leading contributions to the order pa- 
rameter and the magnetic field, recovers the standard GL 
approach, while the next orders in r constitute its exten- 
sion. Such a hierarchy of equations should be solved re- 
cursively, starting from the standard GL equations. We 
derived and studied the equations for the next-to-leading 
contributions to the order parameter A and the mag- 
netic field B. In order to select all relevant terms in 
the case of a nonzero magnetic field, the normal-metal 
temperature Green function was generalized beyond the 
standard Peierls phase approximation to incorporate ad- 
ditional terms up to the order t 2 . The relevant boundary 
conditions were shown to be directly related to the series 
expansion in r for the current. The accuracy of the GL 
theory extended to the next-to-leading order in t was 
tested by comparing the results of the extended formal- 
ism for the uniform order parameter and the critical mag- 
netic field with the corresponding results of the standard 
GL approach and the BCS theory. This demonstrated 
that the validity domain of the GL theory is consider- 
ably increased by the extension. We found very good 
agreement with the full BCS calculations down to tem- 
peratures 0.3-0.4 T c . To illustrate advantages of the con- 
structed extension to the GL formalism, the surface en- 
ergy for the interface between the superconducting and 
normal phase was investigated. We have found, in a semi- 
analytical form, the temperature-dependent correction to 
the value of the GL parameter k at which the surface 
energy becomes zero. Surprisingly, the obtained correc- 
tion is extremely small: it does not exceed 2% even at 
T = 0.3-0.4 T c . This result implies that the boundary 
between type-I and type-II superconducting behavior is 
almost independent of temperature. 

It should be noted that a functional similar to Eq. (|2"7} 
was considered in Ref. [H] in the context of the FFLO 
state (we have an additional term oc 04). However, there 
is a conceptional difference with our work: we focus on a 
series expansion in r and, so, this functional is only the 
initial point to construct such an expansion. Our focus 
on a perturbation theory in r allowed us to make a proper 
selection of all the relevant contributions. It means that 
we did not simply borrow some functional from previous 
papers as the initial step for our study but instead per- 
formed an extensive procedure of microscopic derivations 
(see Appendices) accompanied by an accurate analysis of 
the temperature dependence of each contributing term. 
We proved that the initial free-energy functional given 
by Eq. (|27p contains all the terms that contribute to A 
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and B in the leading and next-to-leading orders in r. 

We also note that though the term oc 0,4 in Eq. (|27j) 
does not produce a pronounced contribution for con- 
ventional bulk superconductors [it is proportional to 
l/(fcpA) 2 , see, e.g., Eq. (|46d|) ]. going beyond the Peierls 
phase approximation can be of importance, e.g., for 
multi-band (subband) materials/systems, where one of 
the relevant bands is characterized by an extremely small 
Fermi momentum but A is determined mostly by other 
bands with large hp (so that the product kpX can be 
even smaller than 1). One of possible examples is single- 
crystalline metallic superconducting nanofilms, where 
different subbands are induced by the quantization of 
the perpendicular electron motion and the energetic po- 
sition of each subband (with respect to the Fermi level) 
can vary significantly with changing nanofilm thickness, 
substrate material etc.- 24 Similar physics can be expected 
for a thin superconducting layer induced by an external 
electric field at the interface between the semiconductor 
(SrTi0 3 and KTa0 3 ) and an electrolyte^ 



Eq. ([AT]l. i.e., 



,1/2 



9 



A d 3 zK a (z), 



7(3) 



9 

r 5/2 

9 



2! 

d 3 zK a (z) { ^^A, 



(A2a) 
(A2b) 
(A2c) 



where A(r) is the scaled order parameter as a function 

of scaled coordinates. Details on calculations of and 

la can be found in textbooks, e.g., in Ref. [IH with the 
result 

/W=a 1 r 1 / 2 A, 1^ = a 2 r 3 / 2 V 2 A, (A3) 

where a\ and ai are given by Eq. (|12l) . To find ii , it is 
first convenient to rearrange Eq. (|A4[) in the form (odd 
powers of V do not contribute) 
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Appendix A: Coefficients in the expansion of the 
gap equation for zero magnetic field 

1. Coefficients at related to the integral kernel 

K a (v,y) 

We start our derivation of the r.h.s. of Eq. (jTTJ) from 
the terms coming from the integral involving the kernel 
K a (r,y) = K a (z) (with z = r - y), i.e., 



I a = - fd 3 zK a (z)A(r + z). 
9 J 



(Al) 



Following the usual practice, this integral is expanded 
in terms of the spatial derivatives of the order param- 
eter A(r), i.e., Eq. ©. Keeping to the mnemonic rule 
formulated in Sec. HH we conclude that working to the 
order r 5 / 2 , it is necessary to incorporate all the spatial 
derivatives up to the fourth order in the gradient expan- 
sion Eq. (J9j> . Due to the symmetry of the kernel K a (z) 
with respect to the transformation z — > — z, the first- and 
third-order derivatives do not contribute. So, we obtain 
only the three relevant terms in Eq. (|lip that come from 



Z (3) = T _ l d s zKi 



'n v n 



o E ^V 2 V 2 J A, (A4) 



with z = {2:1,2:2,2:3} and V = {Vi,V2,V3}. As seen 
from Eq. (|A4[) . two integrals are needed to be calculated, 
i.e., 

J d 3 z K a (z) z 4 and J d 3 z K a (z) z 2 z 2 n (n ^ m). 

Below we assume spherical symmetry, i.e., K a (z) = 
K a (z), with z = \z\. In this case the above integrals 
do not depend on indices n and m. 
As follows from Eq. ([S]), 



7riV(0) 

M r_ y| 



i k,F |r— y | sgncu — i^- 1 r — y | 



(A5) 



with the Fermi wavenumber, N(0) = mkF/('2iT 2 h 2 ) 
the density of states at the Fermi surface and vf the 
Fermi velocity. When inserting Eq. (|A5I) into Eq. ©, 

one can find [^ 0) (r,y) = -^i°2(y,r)] 



K a (z) = gT 



irN(O) 



tFZ 



sinh(2/^ T ) 



(A6) 



with £ T = hv F /(2nT). Then, based on Eq. §M$) and 
replacing T by T c , we obtain 



J d 3 z K a (z) z 4 = - gcrfvp, 
Jd 3 zK a (z) zlzf n = — 



gc h v F (n ^ ffl), 



(A7a) 
(A7b) 
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where c is given by Eq. (fl"2j) . The integrals can easily be 
taken by using (see, e.g., Appendixes in Ref. U8I ) 



J,-i 



dx 



sinh(x^ 



= 2(i-2- £ )r(£)cW (*>i), 



with r(. . .) the Euler gamma-function and £(• . .) the Rie- 
mann zeta-function. Equations (|A4[) and (1A7[) make it 
possible to get 



4 3) = ^ 5/2 ^^ 4 4[E^ + E V^J A, (A8) 



which is reduced to 



/(3) =a 3 r 5 / 2 V 2 (V 2 A), 



(A9) 



orders r 3 / 2 and r 5 / 2 are the following: 



3/2 



r(2) _ r ' 



.9 

5/2 



9 



3 

A |A| 2 /n d3 ^ za, as), (A12b 

»=i 

3 

A yn^^ ^6( z l: z 2,Z 3 ) 



r(3) _ 



x (z 2 - V)A* ((zi+z 3 )- V)A, (A12b) 

r 5/2 , 3 



.9 



A* Jf[d 3 *iKb(*u*2,*a) 

i=X 



r(4) 



r(5) _ 



T 5/2 A 2 



5 2 



x (z x ■ V)A (z 3 • V)A, (A12c) 

3 

/^Qd 3 Zi if h (zi,Z 2 ,Z 3 ) 

^ 1=1 

x(z 2 -V) 2 A*, (A12d) 



t 5 / 2 |A| 2 

<7 2 



3 

/J^[d 3 Zi K b (zi,z 2 ,z 3 



(( Zl • V) 2 + (z 3 • V) 2 )A, (A12e) 



with the coefficient a 3 given by Eq. (IT21 . Now, collecting 
the results for I a 2 ^ and /a 3 '*, we obtain 

I a = aiT 1/2 A + a 2 r 3/2 V 2 A + a 3 r 5/2 V 2 (V 2 A), (A10) 



The contribution given by Eq. (|A12a[) appears already in 
the standard GL domain, and, so, explanations on eval- 
uating I^p can be found in textbooks, see, e.g., Refs. 
and 18. This term reads 



= -6 1 r 3 / 2 |A| 2 A. 



(A13) 



as appears in Eq. (fTTj) . 



2. Coefficients bi related to the integral kernel 

K b (r, Ms) 



Our next step is to calculate the coefficients bi in 
Eq. (|TT|) that are related to the second integral kernel, 
i.e., Kt,(r, {y} 3 ). We start with 



h = - J\\_^ z i K b (z 1 ,z 2 ,z 3 ) 



x A(r + zi)A*(r + z 2 )A(r + z 3 ), (All) 



with Zj = yi - r (i = {1,2,3}) and K b (zi, z 2 , z 3 ) = 
K b (r,{y} 3 ). The integral in Eq. (jAllI) is expanded in 
terms of the spatial derivatives of A(r) and A*(r) by 
using Eq. ©. The terms that contribute to the relevant 



The other terms in Eq. (|A12I) require a more involved 
calculational procedure. Below the details of such a pro- 
cedure are given for 1^ . As to calculations of , 7^ 

(5) 

and 1^ , we restrict ourselves to only basic remarks. 

(2) 

The term 1^ can be written as 



= - T ^T C (1 + 0(t)) E A V„A V to A* 

nrn 

x E / li^^Hi)^^-^) 

u J J=l 

x gW(z 2 - z 3 ) ^I 0) (z 3 ) («i„ + «3»), (A14) 



with ^ 0) (r,r') = g<^ uj (r - r') and Zj = {^1,^2,^3}. 
The integral in Eq. (|A14I) is reduced by invoking the 
Fourier transform and applying the well-known convolu- 
tion theorem provided that we rearrange the polynomial 
in the relevant integrand as 

Z2m (Zln + Z 3n ) = {-Zi m )(-Z ln ) - (-Z lm )z 3n 

-\-[Zlm — Z 2m )( — Zi n ) — \Z\m — Z 2m )z 3n . 



?(0), 
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Then, one finds 
r(2 



jW =t */3T c (i + ( t )) 
:£A V„A V m A*£ / — 

nm " J ^ 

: \ (d m d n — — 



d 3 fc 
(2tt) 3 

1 



( m jfta; — ) ( " ifku + ) ft 2 w 2 



+ 3 



( m iftw + ft- ) ( " iftw + 



iftu; + f fe 7 (iftu; - £ fc ) 2 



with <9 k = ^2, 5 3 }. After straightforward but tedious 
calculations we further obtain 

j(a) = _ r 5/2 Tc (! + 0( r) ) ^ A V„A V„A* 

n 

'd 3 /c ft 4 fc 2 r 2 



E 



(2tt) 3 m 2 L(iftw-f fe ) 4 (^ + &) 2 



(ifuj-^) s (ifuj + ^) 3 y (A16) 

with k = {ki,k2,k^}. Due to the spherical symmetry 
of the term in the parenthesis, the integral in Eq. (|A16p 
does not depend on n so that fc 2 can be replaced by 
fc 2 /3. Then, making use of the standard approximation 

(£ = &) 

iv(o) /de, 



d 3 fc 



(2tt) 3 



one gets 

r(2) 



4ft 2 



r 5 / 2 T c (l + 0{t)) A |VA| 2 ^- JV(0) 



+ 00 



(f + /z)(2g 2 + 2zft^) 
(ft 2 w 2 +£ 2 ) 4 



(A17) 



The terms in the numerator of the integrand proportional 
to an odd power of £ do not contribute. The same is re- 
lated to the terms proportional to uj due to the summa- 
tion over the positive and negative Matsubara frequen- 
cies. So, Eq. (|A17[) is further reduced to 



r(2) 



r 5 / 2 T c (l + G(r))A|VA| 



x 7V(0)^ 



4ft 2 
3m 



T— 



da 



2a 2 



(1 + a 2 ) 4 



(A18) 



where, we recall, huj = 7rT(2n +1). Now, using 



1 



31C(5), / da 



2a 2 



(1 + a 2 ) 4 



with £(. . .) the Riemann zeta-function, we arrive at 



4 2) = -26 2 r 5 / 2 A|VA| 2 , 



with 62 given by Eq. (fl2)) . 



(A19) 



(2) 

Based on the calculations of 1^ , we can further pro- 
ceed with 1^ , 1^ and 1^. The contribution given by 
1^ can be reduced to 

r(3) 



jW _ _ r 5/2 Tc (! + 0( T ))A* (VA) : 



ofe2 1 /■ 1 

^(O)M^y^^ da- =^5, (A20) 

wf 3m ^ |ftw| 5 J (1 + a 2 ) 3 v ' 



>, (A15) with 



+00 



da 



3tt 



(1 + a 2 ) 3 



When making the summation over u, Eq. (|A14[) becomes 
of the form 



4 (3) = -36 2 r 5 / 2 A* (VA) 2 



(A21) 



For ; we obtain 

jW =T 5/2 Tc ( 1 + 0(r) )A 2 v 2 A* 

+ OO 

AT/ „, 2ft 2 v-^ 1 / 3a 2 -1 
^(°)^E^/da— (A22) 



with 



+ OO 



da 



3a 2 - 1 7T 
(1 + a- 2 ) 3 ~ ~8' 



This results in 



At last, 1^ is reduced to 

j(5) =T 5/2 Tc ( 1 + 0(r) )|A| 2 V 2 A 

+ OO 



I, W =-6 2 r 5 / 2 A 2 V 2 A*. 



(A23) 



(1 + a 2 ) 4 



which results in 



il 5 > = -4& 2 t 5 / 2 |A| 2 v 2 A*. 



(A25) 



N ow, ba sed on Eqs. (|ATT]) . (|AT3|) . ([ATS]) . (|MT| . (|Al3)) 
and (|A25[) . we obtain 

4 = - 6!T 3 / 2 |A| 2 A - 6 2 t 5 / 2 [2A |VA| 2 

+ 3A*(VA) 2 + A 2 V 2 A* + 4|A| 2 V 2 A] , (A26) 
with b\ and b 2 defined by Eq. (fT2"l) . 
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3. Coefficient ci coming from K c (r, {y}s) 

The term with the coefficients c\ in Eq. (fTT|) appears 
due to the contribution to A(r)/g given by 

I c = /JJd 3 ^ iC c (zi,z 2j Z3,Z4,z s ) A(r + zi) 

x A*(r + z 2 )A(r + z 3 )A*(r + z 4 )A(r + z 5 ). (A27) 

We need all contributions up to the order t 5 / 2 in Eq. ([TT]) . 
As the leading-order term in the order parameter is pro- 
portional to r 1 / 2 , it is possible to neglect the contribution 
of the spatial derivatives of the order parameter and limit 
ourselves only to the local contribution given by 

/« =r 5 / 2 (l + 0(r)) A|A| 4 

5 

x Jl^zj K c (z 1 ,z 2 ,z 3 ,z 4 ,z 5 ). (A28) 
i=i 

Using Eq. ([5]). performing the Fourier transformation, 
and passing to the integration over the single-electron 
energy, we can find 



j(l) = T 5/2 Tc ^ + A | A |4 N ^ 



y— 



do: 



M 5 J (a 2 + l) 3 ' 



(A29) 



where the integral is equal to 37r/8, see the previous sub- 
section. Evaluating the sum over the Matsubara frequen- 
cies in Eq. (| A29I) . we obtain 



I e = c lT 5/2 A |A| 
where c\ is given by Eq. (IT21 . 



(A30) 



Appendix B: The normal-state Green function in 
the presence of a magnetic field 

As discussed in the main text of the article, construct- 
ing the EGL formalism requires the calculation of the 
normal-state Green function in the presence of a mag- 
netic field with accuracy 0(t 2 ). This means that we are 
not able to rely upon the phase-integral approximation 
of Gor'kov given by Eq. (|18p . where the classical particle 
trajectory is assumed to be a straight line. To go beyond 
this approximation, we need to take r-dependent devi- 
ations from such a linear trajectory. A natural way of 
doing so is based on the following representation for the 
single-particle propagator (see, e.g., Ref. [l9t): 



0(°)(rf,r't') = F(t-t') e** 
where S c \ is the classical action 



S cl = J'^f ds + - c J A(q)-dq, 



(Bl) 



(B2) 



where the integrals are taken along the classical trajec- 
tory C that satisfies the equation of motion 



m c 



(qxB) 



(B3) 



with the boundary conditions q(t') = r', q(f) = r. We 
are interested in the systematic corrections to Gor'kovs 
cikonal approximation and, so, it is of convenience to 
recast S c \ in the form: 



Sr] — St 



Gor 



■ — r' \ 2 



t-t' 



ds- 



B(q)-dS, 

(B4) 

where Scor is the classical action along the straight line 
connecting r' and r, i.e., 



S'cor — 



m(r — r') 1 
2{t-t') 



dq, 



(B5) 



which is the basis of the Gor'kovs approximation for the 
normal-state Green function in a magnetic field. The in- 
tegral in Eq. (|B2I) is over the surface S that is bound by 
the loop (9£ consisting of two parts, i.e., the classical tra- 
jectory C from r' to r and the straight line connecting r 
and r'. The orientation of the <9£ is positive with respect 
to the surface. 

To simplify the further analysis, we introduce the de- 
composition q = qx + q 2 and recast Eq. (jB3[) as follows: 



qi 



m c 

e 

m c 



qi x 



B(r)), 



(B6a) 



^ = — fa x t 5 ^) - B ( r )]) + — fa x B ( £ l))' 

(B6b) 



with the boundary conditions qi(£') = r', qx(t) = r and 
q2(i) = q2(i') = 0. This decomposition is such that 
qi is the solution for the uniform (not dependent on q) 
magnetic field B(r). This rearrangement is convenient 
because, as seen below, the spatial derivatives of B(r) 
contribute to the propagator only in the order r 5 / 2 . Now, 
we set B(r) = B(r)e z , with e z the unit vector in the z- 
direction. The corresponding r-expansion for qi can be 
found from straightforward calculations with the result 
given by (for the moment, we put r' = 0, t' = 0) 

q lx = x s/t + t yQ(p(s) + t 2 xn 2 X (s) + 0{t 3 ), (B7a) 
q ly = y s/t - t xti^{s) + t 2 yti 2 X (s) + <D{t 3 ), (B7b) 



(Liz 



z s/t, 



(B7c) 



where Cl = \e\B(r)/mc, with the scaled magnetic field 
B(r) — -B(r) (in other words, f2 stands for the scaled 
cyclotron frequency). In addition, qi = {q\ x , qi y , qi z } 
and r = {x, y, z} in Eq. (|B7|) . and (p(s) and x( s ) are 
given by 



<p(s) 



1 



X(s) 



1 



t 

3s 



2s- 
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where, as seen, ip(t) — x(t) = 0. It is worth noting that 
B(r) does not depend on r, see Eq. (|2"9")l . 

The solution of Eq. (|B6b|) is more complicated because 
it involves the spatial derivatives of B(r). Let us first con- 
sider a simplified variant when the magnetic field is gen- 
erally parallel to the z-axis while its absolute value varies 
with position, i.e., B(q) = B(q)e z . Then, Eq. (|B6b[) can 
be rearranged as 

q 2 = — [B^) - B(r)] (q< 0) x e.) + 0(r 5 / 2 ), (B8) 
m c 

with q^ =lim T _ ) .oqi) i-e., q^ = {xs/t,ys/t, zs/t}. We 
can further simplify this equation by making use of 

B(ctf ) )-B(r) = 

= r 3 / 2 i (r • V)B(r) + r 2 ^(r ■ V) 2 B{v) + G(r 5 / 2 ), 

where scaled derivatives are introduced, i.e., V = 
T — 1 / 2 V, and the operator V acts only on the magnetic 
field. From the resulting equation, we obtain 

V 



to = |S( S ,r)+G(r 5 / 2 ), 

5/2 1 



Q2y 



t 



S( S ,r) + G(r 5 / 2 ), 



(B9a) 
(B9b) 

q 2z = 0, (B9c) 
where S(s, r) satisfies the differential equation 



E(s, r) = J- [r 3 / 2 f (r ■ V)B(r) + r 2 ^(r ■ V) 2 B(r)] , 



supplemented by the boundary conditions S(0, r) = 
S(i,r)=0. 

Now, we have everything at our disposal to calculate 
Sci — Sqoi and, so, to analyze the corrections to the 
approximation employed by Gor'kov. Using Eqs. (|B2I) . 
(Ell), (lB7l> and (IBSl). we find 



Sol - S Gor = --r 2 (x 2 + y 2 )tt 2 t + G(r 5 / 2 ). (BIO) 

It is remarkable that the spatial derivatives of the mag- 
netic field contribute to S c i — Sq oi: only in orders higher 
than r 2 . In particular, when considering the surface inte- 
gral in Eq. (|B4[) , Bar and the contribution of the spa- 
tial derivatives of the magnetic field to the surface inte- 
gral, taken in its lowest order, is proportional to t 3 / 2 . So, 
the resulting product is of the order r 5 / 2 , which means 
that it does not make a contribution to the first term in 
the right- hand- side of Eq. (|B10|) . In turn, calculating the 
kinetic energy we find 



qi 



5( S ,r) 



t 



(yqix - xq ly ^j + 0(t 5/2 ), 



which, taken together with q\ x — x/t and q\ y = y/t, 
makes it possible to conclude that the spatial derivatives 
of the magnetic field can contribute to S c \ only to the 
order t 5 / 2 . 



Let us make a few remarks about Eq. (|B10[) . Though 
this result was derived under the simplified assumption 
B(q) = B(q)e z , it is general and holds even in the pres- 
ence of spatial variations in the direction of the magnetic 
field. Thiscan be seen from the following. Based on our 
above consideration, we expect that the two lowest orders 
contributing to B(q^°' ) ) — B(r) are t 3 / 2 and r 2 . Then, 
Eq. (|B6b|) can be reduced to 

q 2 = — (qi 0) x [B(ql O) )-B(r)]) + 0(T 5 / 2 ), (Bn) 
mt 

whose solution reads (i — {x, y, z}) 



met 



^e ijkrj r k (s,r)+0(T 5 / 2 ), (B12) 



with Ejmk the permutation tensor and Bj(q^) — Bj(r) = 
Ti, where Tj is taken in the two lowest orders in r (with 
the boundary conditions Yj(0,r) = Ti(t, r) = 0). What- 
ever Ti, it does not make a contribution neither to the 
surface integral nor to the kinetic term in each order lower 
than r 5 / 2 . In particular, in the kinetic term we obtain 



ijk 



where the second term is simply equal to zero. We also 
remark that detailed calculations make it possible to find 
that 



r 3/2 



q2 



6m c 



1-- (r.V)(B(r)xr)+G(r 2 ), (B13) 



which means that the above assumption about the two 
lowest contributing orders in B(q^°' ) ) — B(r) is fully cor- 
rect. 

The only thing remaining is to specify the quantum- 
fluctuation factor F(t) (t' = is still of use). In the 

3 /2 

Gor'kov approximation F(t) = ( 2 ^iht ) but * s n0 ^ 
longer the case in the order r 2 and higher. Within accu- 
racy 0(t 2 ) we obtain 

. , / m \ 3 / 2 / T 2 Cl 2 ,\ 

which is nothing else but the fluctuation factor for the 
propagator in a uniform magnetic field expanded in r 
(see, e.g., Ref.EI). 

Now, based on Eqs. jB4|, (|B5j). (|BT0]> and (|BT4|) and 
making the usual imaginary-time substitution t — > —it, 
we arrive at the following expression for the temperature 
single-particle Green function in the presence of a mag- 
netic field (t — > t — t' and r —> r — r'): 



24 



(t-f) 



>\2 



+ ^(r-r')l(<-<')]+0(r 5/2 )}, (B15) 
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where (r — r')j_ is the component of the vector r — r' 

perpendicular to B(r) and <?Qo r (rt, r't') is the Gor'kov 
approximation for the normal-state Green function given 
by Eq. ([TBI) . Now, switching to the Fourier transform of 
the Green function given by Eq. (IB15|) we find Eq. ([20]) . 
As seen from Eqs. (l20l) and (|B15|) . the spatial deriva- 
tives contribute to the r-expansion of the temperature 
normal-state Green function only in the order r 5 / 2 and 
higher. However, to construct the EGL formalism (up to 
the order r 3 / 2 in the order parameter), we need to know 
\r, r') only up to the order r 2 . 



orders in r, we obtain the following contributions: 



Appendix C: Expansion for the gap equation in the 
presence of a nonzero magnetic field 



In Appendix [X] we gave the details of the calculations 
for the coefficients appearing in the r-expansion of the 
gap equation. It is of great importance to specify com- 
plications that appear when generalizing the procedure to 
the case of a nonzero magnetic field, i.e., for the normal- 
state temperature Green function given by Eq. ([20)) . 



1. Terms related to the integral kernel K a (r,y) 

When switching to a nonzero magnetic field, Eq. (|A1[) 
in Appendix [A] can be rewritten in the form 

I a =- /dV^(r,y)A(y) 
9 J 

= lim - ( d 3 z Q Jr, z)A(r + z, r'), (CI) 

r'^r g J 



where the "two-point" order parameter A(r, r') is defined 
by Eq. (HQ) and 



e 2 B 2 (r) 



) a (r, z) = K^ B =o(z) - gT c 



24m 2 c 



(C2) 



with K a B=o( z ) the zero- magnetic field kernel K a given 
by Eqs. © and (JSJ). Then, introducing the expansion of 
A(y, r') in powers of z = y — r and collecting the relevant 



± a 

r(2) 

± a 

r(3) = 

a 

/i 4 > = 



T V2 _ /■ 

lim A(r,r') d 3 z K a . B=0 {z), (C3a) 

lim d 3 zK a , B=0 (z) { o , rJ A(r,r'), 



5 >' r 

-5/2 

lim /d 3 z^ aiB=0 (z) 

.9 r '^ r 



2! 

(z-V r ) 4 



4! 



(C3b) 

A(f,f), 

(G3c) 



12m z C z r'->r 



d"* £<Uo(-)(# - ^mzi^)^°L (z) 



(C3d) 



After integrating over z (for liVfU 33 see details in 
Appendix |A1 for I a see the discussion below), Eqs. (IC3|) 
are reduced to 

/(!) = air 1 / 2 lim A(f,r'), (C4a) 

r'— >r 

/( 2 ) = a 2 r 3 / 2 lim V 2 A(f , r'), (C4b) 

r'— >r 

= a 3 r 5/2 lim V 2 (V 2 A(r, f')) , (C4c) 

r'— «• 

4 4) = -a 4 r 5 / 2 B 2 (r) lim A(f,f'), (C4d) 

r'— »r 

where the coefficients a±, a?, a 3 are given by Eq. (j!2[) and 
04 is defined in Eq. (|251) . 

As already mentioned in Appendix El the first two 
terms, i.e., Ia and I a 2 ^ appear even in the standard GL 
domain and, so, the details of calculating these contribu- 
tions are well-known from textbooks; 2 *^ The results are 
given by Eqs. (|22aj) and (|22b|) . and the only difference 
from the standard GL theory is that the coefficients a\ 
and 02 contains now extra terms of the order t 2 . Calcu- 

( 3) 

lating I a is a more involved and complicated task and, 

(3) 

so, the basic details of calculating I a are outlined below. 
We remark that 



lim (V 2 ) 2 A(r,r') 



= lim 

r'— H 

with <&' r (v,v') = V r $(r,f'), where 



V r -— $;(r,r')) ] A(r), (C5) 



$(f,f')= / A dq, 



(C6) 



here the integral is taken along the straight line connect- 
ing the points r' and r. Expanding A(q) in powers of 
q — r' we can rewrite $(r, r') in the form 



^ (n + 1)! 

n— 



a— r — r' 



(C7) 
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As seen, lim V r $(r,r') = A(r) and, so, one could ex- "two-point" order parameter, i.e., given by Eq. (|2~Tj) . as 



pect that lim (V?) = (D 2 ) for n = 1,2 in Eqs. (ICibl) 

r'— sr 

and (|C4c[) . However, this is true only for n = 1 (i.e., for 

la ) and does not hold for n = 2: the limiting proce- 
dure and differentiating do not commute in general. In 
particular, straightforward but tedious calculations show 
that 



lim (V 2 ) 2 (A(f) e-ff^') 

r'— >r \ 



D 2 (D 2 A) - ViA(ViV,-A J - - V 2 A 4 ) 



37ic 



4e 
ft 2 



V,.<1 ; (V,.l, V,.l ; ; 



--^.(V^-V.-V^)], (C8) 

with A = {A 1 ,A 2 ,A 3 } and V = {Vi,V 2 ,V 3 }. The 
above expression results immediately in Eq. (I22c[) , when 
keeping in mind that VjA; — \?jAi = J2k £i J k -Bk, with 
Sijk the permutation tensor and B = {i?i, i3 2 , -B 3 }. 

Concluding this subsection, we note that integration 
over z in 1^ (together with the accompanying summa- 
tion over the Matsubara frequencies) is similar to that for 
Ia^ with i = 1, 2, 3. Invoking the Fourier transformation 
and converting z 2 j_ into the corresponding derivatives of 
the Green functions with respect to the single-particle 
energy, for ilf 1 we find 

7 W = - T ^T c (l + 0(r))g^li m A(f ) f') 



yI—{- 

^ J (27t) 3 I (i 



12m 2 C 2 r'-s-r 



(2tt)3 l(ihLj + Z k )(ihu-£ k )3 



with dk = {di, 82,83}. This is reduced to 
4 4) = -r 5 / 2 T c (l + OM) lim A(r,f') 



h 2 e 2 B 2 (r) 
9m 2 c 2 



(a 2 + l) 3 ' 



Here the integral is equal to 7r/8, and for the sum 
over the Matsubara frequencies we obtain J2 U VI^M 3 = 
7C(3)/(47r 3 T3)(l + 0(T)). This gives Eq. (Icfdl . 



2. Terms related to Kb(r, {y}3) 

The second generation of the terms appearing in the 
r.h.s. of the field-modified Eq. (J7J comes from Eq. (|A11|> . 
which can now be rewritten in terms of the auxiliary 



I b = lim - f]d 3 Zi Q b (r, zi,z 2 ,z 3 ) 

x A(r + z 1 ,r')A*(r + z 2 ,r')A(r + z 3 ,r'), (Cll) 

where 

Q b (r, zi,z 2 ,z 3 ) = isr 6)B= o(zi,z 2 ,z 3 ) e s (G12) 
and £ = £(r, Zi, z 2 , z 3 ) is given by 
E(r,zi,z 2 ,z 3 ) = 

= $(r,r + z 2 ) + $(r + z 2 ,r + zi) + $(r + zi,r) 

+ $(r, r + z 2 ) + $(r + z 2 , r + z 3 ) + $(r + z 3 , r), 

(C13) 



with $(r,r') given by Eq. (|C6]l . Equation (|C13|) can be 
rewritten in the form 

£ = yy"B(r)-dS!+ JjB(r)-dS 2 +0(T 3 / 2 ), (C14) 

where the boundary of the surface Si consists of the three 
line segments that connect the points r, r+zi, and r+z 2 , 
i.e., r — » r + zi — » r + z 2 — > r; and S 2 is bounded by the 
three line segments between the points r, r + z 3 , and 
r + z 2 , i.e., r —> r + z 3 — >• r + z 2 —> r. Linear boundaries 
make it possible to analytically calculate £(r, Zi, z 2 , z 3 ), 
i.e., 

E = r(B(f) • [(zi x z 2 ) + (z 3 x z 2 )]) + 0(r 3 / 2 ). 

(C15) 

Now, we have everything at our disposal to collect the 
relevant terms up to 0(t 5 / 2 ) in Eq. (fTTj) . Based on the 
consideration in Appendix [A~l one obtains 



7 b (1) = -&ir 3 / 2 lim |A(r,r')| 2 A(f,f'), 



(C16) 



4 (2) = -b 2 r 5/2 lim 



2A(r,r') |V r A(r , f')| 2 

+ 3A*(f,f')(V r A(f,f')) 2 +A 2 (r,r') 
x V 2 A*(r, f') + 4| A(r , f')| 2 V 2 A(r, f') 



(C17) 



which results in Eqs. (|23| (here V can safely be replaced 
by D in the relevant expressions). 

It may seem that there is one more term of the or- 
der t 5 / 2 coming from the kernel Kf, in the presence of a 
nonzero magnetic field, i.e., 



r(3) _ T 
1 b ~ 



5/2 



3 

A| 2 A Hd 3 z,^ fc(1 (z 1 ,z 2 ,Z3) 

x (fi(r). [(zixz 2 ) + (z 3 xz 2 )]). (C18) 
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However, we immediately find that I b — due to 
K b ,B=o(zi, 22,23) = K b)B =o(-zi, -z 2 ,-z 3 ). 

Concluding Appendix [Cj we note that the only term 
appearing in Eq. ([lip from the integral with the kernel 



K c is proportional to |A| 4 A and, so, does not change its 
form in the presence of a nonzero magnetic field (here 
corrections from B^O can appear only in higher orders 
in t). 
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